#include <iostream>
#include <fstream>
#include <stdexcept>

#include "om_matrix.h"

#include "minimization.h"
#include "operators.h"

class FunctionOfInterest
{
public:
	typedef om::vector<3> vector;
	typedef om::matrix<3> matrix;

	FunctionOfInterest(const double alpha) : alpha_(alpha) {}

	double operator() (const vector& x) const
	{
		return ( pow(x[0]-x[1], 2.0) + pow(x[0], 4.0)*pow(x[1]-0.5, 2.0) + 2.0*pow(x[2]-x[1]+alpha_, 2.0) );
	}

	vector gradient(const vector& x) const
	{
		vector g;
		g[0]=4.0*pow(x[0], 3.0)*pow((x[1]-0.5), 2.0) + 2.0*(x[0]-x[1]);
		g[1]=(0.0-4.0)*(x[2]-x[1]+alpha_) + 2*pow(x[0], 4.0)*(x[1]-0.5) - 2.0*(x[0]-x[1]);
		g[2]=4.0*(x[2]-x[1]+alpha_);
		return g;
	}

	matrix hessian(const vector& x) const
	{
		matrix h;
		h[0][0]=12.0 * pow(x[0], 2.0) * pow((x[1]-0.5), 2.0) + 2.0;
		h[0][1]=8.0 * pow(x[0], 3.0) * (x[1]-0.5) - 2.0;
		h[0][2]=0.0;
		h[1][0]=8.0 * pow(x[0], 3.0) * (x[1]-0.5) - 2.0;
		h[1][1]=2.0 * pow(x[0], 4.0) + 6.0;
		h[1][2]=-4.0;
		h[2][0]=0.0;
		h[2][1]=-4.0;
		h[2][2]=4.0;
		return h;
	}

private:
	double alpha_;
};

struct Configuration
{
	double alpha;
	FunctionOfInterest::vector start;

	Configuration(int argc, const char** argv) : alpha(0.0)
	{
		if(argc<2) { throw std::invalid_argument("Missing arguments"); }
		std::ifstream input(argv[1], std::ios::in);
		if(!input.good()) { throw std::runtime_error("Failed to open input file"); }
		input >> alpha;
		start=FunctionOfInterest::vector::from_stream(input);
	}
};

int main(int argc, const char** argv)
{
	try
	{
		Configuration conf(argc, argv);
		const FunctionOfInterest function(conf.alpha);
		ConvergenceChecker<FunctionOfInterest> checker(om::epsilon(), 50);
		const NewtonianDirector<FunctionOfInterest> director;
		const GoldenSearcher<FunctionOfInterest> stepper(1.0, om::epsilon());
		minimization<FunctionOfInterest>::minimize(function, conf.start, checker, director, stepper);
	}
	catch(const std::exception& e)
	{
		std::cerr << "Exception: " << e.what() << "\n";
	}

	return 0;
}
